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Abstract 

Background: Cryptic population structure can be an indicator of incipient speciation or historical processes. We 
investigated a previously documented deep break in the mitochondrial haplotypes of Heliconius eroto chestertonii 
to explore the possibility of cryptic speciation, and also the possible presence of endosymbiont bacteria that might 
drive mitochondrial population structure. 

Results: Among a sample of 315 individuals from 16 populations of western Colombia, two principal mtDNA 
clades were detected with 2.15% divergence and we confirmed this structure was weakly associated with 
geography. The first mtDNA clade included 87% of individuals from northern populations and was the sister group 
of H. eroto members of Andes western, while the second clade contained most individuals from southern 
populations (78%), which shared haplotypes with an Ecuadorian race of H. eroto. In contrast, analysis using AFLP 
markers showed H. e. chestertonii to be a genetically homogeneous species with no association between 
mitochondrial divergence and AFLP structure. The lack of congruence between molecular markers suggests that 
cryptic speciation is not a plausible explanation for the deep mitochondrial divergence in H. e chestertonii. We also 
carried out the first tests for the presence of endosymbiontic bacteria in Heliconius, and identified two distinct 
lineages of Wolbochio within H. e. chestertonii. However, neither of the principal mitochondrial clades of/-/, e. 
chestertonii was directly associated with the patterns of infection. 

Conclusions: We conclude that historical demographic processes are the most likely explanation for the high 
mitochondrial differentiation in H. e. chestertonii, perhaps due to gene flow between Cauca valley H. e. chestertonii 
and west Pacific slope populations of H. eroto. 



Background 

Sequences derived from the mitochondrial genome are 
commonly used both in species delimitation and histori- 
cal phylogeography. For example, deep divergence in 
mitochondrial DNA sequences (mtDNA) between 
related individuals is often taken as evidence for the 
existence of cryptic species [1-4]. The discovery of cryp- 
tic species-level variation has important implications for 
characterising biodiversity and for studies of speciation. 
Nonetheless, it is now well recognised that inference of 
evolutionary processes and species boundaries from 
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mitochondrial sequences alone can be problematic [5,6]. 
For example, similar patterns of divergence can be due 
to host-parasite interactions, whereby selection leads dif- 
ferent molecular markers to show different histories [7]. 
Mitochondrial lineages may also be retained through 
admixture between divergent species or populations [8], 
or perhaps due to unusual population structures [9,10]. 
Indeed, in diverse tropical radiations, mtDNA 'barcod- 
ing' may perform rather poorly as a species identifica- 
tion tool [6]. It is therefore of general interest to pursue 
individual cases of deep mtDNA divergence in order to 
determine how often such divergence is indeed an indi- 
cator of cryptic species-level variation. 

Heliconius butterflies are an excellent ecological and 
genetic system for studying speciation. These unpalata- 
ble butterflies are recognized for their diversity of wing 
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color patterns associated with mimicry [11,12]. The clas- 
sical example of this adaptive radiation occurs between 
the comimetic species H. erato and H. melpomene. 
These butterflies co-occur in Central and South Amer- 
ica and show convergent changes in their color pattern. 
They are represented by more than 20 different geo- 
graphic forms which are considered subspecies [13]. 
Neutral molecular markers show geographic structure 
among subspecies of H. erato and H. melpomene, and 
most named races fall within a particular geographic 
clade [13-16]. However, in H. erato two forms: H. e. 
hydara and H. e. chestertonii are polyphyletic in the 
mitochondrial phylogeny [13,16]. Individuals of H. e. 
chestertonii fall into two distinct mtDNA clades that 
show over 2% divergence, with no clear biogeographic 
explanation [17]. H. e. chestertonii is found in the wes- 
tern Colombian Andes on disturbed, dry habitats and 
forms a hybrid zone with the geographically closest sub- 
species: H. e. venus. Although H. e. chestertonii is a 
member of the erato clade, it has an unusual wing color 
pattern compared with the characteristic red/yellow/ 
black pattern of H. erato. H. e. chestertonii has an irides- 
cent and melanic forewing while the hindwing displays a 
broad yellow band. The co-mimic for this wing color 
pattern is H. cydno weymeri f. gustavi, a member of the 
H. melpomene clade. 

Rapid evolution in the early stages of the speciation 
process is expected to lead to incongruence between 
morphological and molecular markers [14,18,19]. How- 
ever, in H. e. chestertonii the deep divergence in the 
mitochondrial haplotypes is not readily explained 
through ancestral polymorphism or hybridization. Such 
mtDNA division has not been observed in any other 
Heliconius taxon to date. One possible explanation 
might be that structure is due to the presence of endo- 
symbiotic bacteria [19-21]. Wolbachia are intracellular 
bacteria and infect numerous species of arthropods and 
nematodes [22]. Interactions between these microorgan- 
isms and their eukaryotic hosts often has consequences 
for host reproduction, leading in some cases to breaks 
between populations or species [23,24]. Wolbachia are 
inherited maternally, so their evolutionary fate is tightly 
linked to that of the mitochondrion [25]. Furthermore, 
hybridization between H. e. chestertonii and its nearby 
relative, H. e. venus, produces partially infertile eggs 
[26]. Hybrid sterility of this form can also be generated 
by endosymbionts where parental populations are 
infected by different strains. 

However, an alternative hypothesis for the mtDNA 
break within the continuous distribution of H. e. chester- 
tonii is cryptic reproductive isolation, unrelated to endo- 
symbionts. Given the emphasis on wing color pattern in 
Heliconius speciation, a potential case of cryptic specia- 
tion would be of considerable interest. To investigate 



these possibilities we first extend previous mtDNA sam- 
pling to better document the distribution of mtDNA 
lineages. In order to test for cryptic speciation, we have 
then complemented these data with nuclear Amplified 
Fragment Length Polymorphisms (AFLP) markers, to 
provide a comparison with nuclear biparentally inherited 
markers. Finally, we have tested for a variety of endo- 
symbiotic bacteria to investigate whether mtDNA struc- 
ture could be a result of patterns of infection among 
populations. 

Results 

Mitochondrial relationships 

Our expanded sampling of mitochondrial sequences con- 
firms the deep divergence within H. e. chestertonii, with a 
similar topology obtained with either Bayesian or parsi- 
mony analyses (Figure 1). As expected, both taxa studied 
here are members of the El clade, containing H. erato 
species from Costa Rica, Panama, Colombia and west 
Ecuador [16]. H. e. venus forms a monophyletic group 
together with some members of H. e. hydara, while indi- 
viduals of H. e. chestertonii showed a complex relation- 
ship with other subspecies and is a polyphyletic taxon. 
Three principal clades are observed. The first primarily 
includes individuals collected in populations north of the 
hybrid zone in Calima River Valley (Figure 1). These hap- 
lotypes form a sister group to the other subspecies of the 
western clade El: H. e. petiverana, H. e venus and H. e. 
hydara (Figure 1). The second principal clade includes 
the west Ecuadorean race H. e. cyrbia and primarily indi- 
viduals of H. e. chestertonii collected in southern popula- 
tions (Figure 1). The third clade, with only three 
northern individuals, is sister to all races of H. erato in 
the western El clade and south clade of H. e. chestertonii. 
The different clades involving members of H. e. chesterto- 
nii has strong bootstrap support (100%) and a high pos- 
terior probability (1) (Figure 1). A haplotype network 
analysis revealed the same relationships among popula- 
tions (Additional file 1). Three main star-like species/ 
geographical clusters, H. e. venus, H. e. chestertonii-North 
and H. e. chestertonii-South. A topology test showed that 
monophyly of H. e. chestertonii could be strongly rejected 
(P < 0.0001 Shimodaira-Hasegawa test). 

Genetic analysis of mitochondrial DNA (mtDNA) within 
populations, showed that H. e. chestertonii is more poly- 
morphic than H. e. venus (0 W - 0.01 and 6 W - 0.003 
respectively, Table 1). We then subdivided H. e. chesterto- 
nii to compare polymorphism within (1) northern and 
southern populations and (2) the two principal clades 
obtained in the phylogenetic tree. The northern and 
southern populations had similar levels of DNA poly- 
morphism (0 W = 0.01 and 0.009 respectively) with no 
fixed difference between groups and 35 shared mutations 
(Table 1). Nonetheless, the comparison between clades 
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Figure 1 Sampling sites of Heliconius erato chestertonii and H. e. venus and bayesian tree for the Clade El of Heliconius erato. Symbols 
correspond to population of origin for the individuals of H. e. chestertonii, south (gray) or north (black) and H. e. venus (white). The tree is based 
on mitochondrial genes of Cytochrome oxidase subunits I and II, leucine-tRNA and has the same topology as the Parsimony analysis. Posterior 
probability and bootstrap support are indicated on branches of the principal clusters. Individuals of the different subspecies are identified by bars 
(right). 
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Table 1 Description of mtDNA polymorphism among 
populations of H. e. chestertonii 



Species/Population/Clade 
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Atunrela (AT) 


27 


5 


0.01 196 


13.73 


0.0103 


Miravalle (MV) 


33 


6 


0.01045 


12 


0.01259 


Montahitas (MO) 


4 


4 


0.00151 


1.73 


0.00153 


La Cumbre (CU) 


33 


6 


0.01341 


15.4 


0.01259 


Saladito (SA) 


24 


4 


0.01238 


14.2 


001004 


Pance (PA) 


24 


3 


0.00697 


8 


0.00916 


Villa Colombia (VC) 


3 


3 


0.00145 


1.66 


0.00114 


Calima River valley (CV) 


29 


5 


0.01096 


12.59 


0.00822 


H. e. chestertonii (CV) 


26 


2 


0.00647 


7.42 


0.00924 


H. e. venus (CV) 


20 


3 


0.00498 


5.71 


0.00711 


North populations 


52 


21 


0.01206 


13.829 


0.01017 


South populations 


42 


16 


0.0077 


8.83 


0.00889 


North Clade 


9 


9 


0.00125 


1.43 


0.00189 


South Clade 


23 


18 


0.00279 


3.19 


0.00456 



The symbols in the table represent: number of segregating sites (S), number 
of haplotypes {h), nucleotide diversity per site (n), average number of 
differences between pair of sequences (/c) and genetic diversity per site {0w). 



showed 17 fixed differences and just three shared muta- 
tions (Table 1). The Analysis of Molecular Variance 
(AMOVA) showed that the maximum variability is within 
populations (63.65%) and only 18.22% of the total variabil- 
ity was due to geography, and a further 18.13% explained 
by species (groups in the AMOVA). When the species 
were analyzed separately, we found higher genetic struc- 
ture among localities of H. e. chestertonii than between H. 
e. venus populations (F ST = 0.2775, P < 0.01; F ST = 0.1197, 
P > 0.5; Additional file 2). 

AFLP Analysis 

Of the 327 loci examined, we found that 318 were poly- 
morphic, representing 97% of loci. Genetic diversity sta- 
tistics estimated from AFLP data showed that 
populations are differentiated within and among species 
(Additional file 3, P < 0.001). Both H. e. venus and H. e. 
chestertonii had values of genetic diversity of 27.6%, 
while the hybrid zone populations (Calima River Valley) 
were more diverse than in any allopatric locality (H. e. 
chestertonii (CV) = 29% and H. e. venus (CV) = 30%). 

The AFLP-based bayesian population assignment test 
identified two clusters which correspond to H. e. ches- 
tertonii and H. e. venus (Figure 2a, L = -20596.26). We 
compared the mean likelihood values for all runs with 
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Figure 2 Structure Analysis, a. Graphical representation of results 
obtained from Structure. Upper black and gray bars represent the 
phenotype of each individual, H. e. chestertonii and H. e. venus 
respectively. The colours represent the Bayesian clusters when the 
analysis was carried out with K = 2 (upper, L = -20596.26) and K = 3 
(lower, L = -20852.973) and correspond to H. e. venus (yelow) and H. 
e. chestertonii (violet and pink). The lower bars and letters show the 
population origin of individuals (for description of locations see 
Table 3). b. The results for K = 3 including only the individuals used 
in mtDNA analysis {L = -12353). Lower black and gray bars represent 
individuals in the southern and northern clades. 

V . J 

the statistic AK and confirmed that K = 2 is the best 
estimate for this data. Additionally, the PCA analysis 
gave a consistent result, with 56% of the variance in the 
data explained by differences between species (Addi- 
tional file 4). In addition, all the populations used in this 
study showed significant genetic structure for AFLP 
markers (Figure 3 and Additional file 2). The AMOVA 
showed that 48% of the total variation in the AFLPs was 
explained within populations and the 42% among 
groups, which correspond to species in our analysis. 
When we excluded H. e. venus and repeated the 
AMOVA grouping data for northern and southern 
populations of H. e. chestertonii, only 9% of the variation 
was accounted for by geography. In order to test for 
congruent differentiation between AFLP and mtDNA 
markers, we ran a Structure analysis using only those 
individuals with data for mtDNA and assuming K = 3 
(putatively representing H. e. venus and the two H. e. 
chestertonii clades). This analysis showed no evidence 
for any population structure associated with the mtDNA 
clades (Figure 2b). 

Interactions with endosymbionts 

PCR analysis of 307 wild caught H. e. chestertonii did 
not identify any positive results for Ricketsia and Spiro- 
plasma infection. All screens were run with a positive 
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Figure 3 Matrix of pairwise F ST . Upper and lower matrices show 
the F ST values for mtDNA and AFLP markers respectively, between 
each population. H. e. chestertonii (black) and H. e. venus (red) 
localities are provided in Table 3. In the Calima River Valley hybrid 
zone, individuals were separated by phenotype: H. e. chestertonii 
(CV) or H. e. venus (CV*). 



control of DNA from known infected individuals of 
other insects, so this failure is unlikely to be due to 
technical PCR error, although we cannot completely 
rule out the presence of a divergent strain of one of 
these bacterial taxa. In contrast, 7% of individuals 
screened (10 females/ llmales) were infected with Wol- 
bachia. There are two major clades within Wolbachia, 
both of which are known to infect insects. Interestingly, 
the sequences derived from our samples fell into both 
clades, with 19 in clade-B and the remaining two in 



clade-A (Table 2). The presence of this endosymbiont is 
almost exclusive to H. e. chestertonii, with just one 
female H. e. venus from the hybrid zone infected. Most 
of these butterflies were collected in northern popula- 
tions of H. e. chestertonii (85%). However not all indivi- 
duals of the same population were infected and some 
localities did not show the presence of Wolbachia. 

Discussion 

It is generally accepted that genetic differentiation 
between subpopulations can lead to the formation of 
new reproductively isolated species over time [27]. Phy- 
logeographic analysis can be useful in identifying both 
cryptic species and incipient subpopulations on the way 
to becoming new species [2,3]. The analysis of this cryp- 
tic population structure can show how genetic, beha- 
vioral and ecological processes have acted during the 
earliest stages of speciation [27]. Studies that include 
independent sources of evidence, such as morphological 
comparisons, reproductive biology and phylogenetics, 
are necessary to understand the history of diverging 
lineages and to resolve species identification [28]. 

In this study, we first confirmed the unusual degree of 
mitochondrial divergence within H. e. chestertonii. Our 
analysis was based on a broad sampling of H. e. chester- 
tonii which enabled the confirmation of two principal 
groups of mitochondrial haplotypes (Figure 1 and Addi- 
tional file 1). Our broader sampling showed that these 
clades are not completely associated with the geographic 
distribution of the populations, as had been suspected 
previously (Figure 2b) [17]. The principal haplotype 



Table 2 Distribution of Wolbachia within populations of 
H. e. chestertonii and H. e. venus 



Population 
(Individuals 
analized) 


Negative 
diagnostic 
(Females/Males) 


Clade-A 
(Females/ 
Males) 


Clade-B 
(Females/ 
Males) 


Marsella (33) 


9/18 




3/3 


Trujillo (13) 


2/10 




0/1 


Yotoco (4) 


1/3 






Buenos Aires (22) 


6/10 


0/1 


1/0 


Carbonero (14) 


6/6 




1/1 


Caimital (15) 


0/9 




3/3 


Atucela (20) 


4/16 






Miravalle (21) 


8/12 




0/1 


Montahitas (20) 


10/9 




0/1 


La cumbre (14) 


4/10 






Saladito (10) 


0/9 




1/0 


Pance (28) 


11/17 






Villa Colombia (19) 


6/13 






Calima River valley 
(44) 


12/31 


170 




Juanchaco (29) 


6/23 






The asterisk indicates the only infected individual of H. e. venus. 
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groups (Figure 1 and Additional file 1) were on average 
2.15% divergent, similar to that estimated in previous 
studies [13,15,17]. This is considerably more divergent 
than within any other race of H erato, apart from the 
trans-Andean biogeographic break which occurs within 
the distribution of the race H e. hydara as described 
above. However, here there is no clear biogeographic 
context for this population structure. Nonetheless, the 
haplotypes do show some geographic structure, with 
87% of northern individuals being Clade 1 and 78% of 
southern individuals Clade 2. 

Cryptic speciation 

If cryptic speciation was occurring, whereby two species 
are considered as one based on wing pattern morphol- 
ogy, we would expect the mtDNA haplotypes to be 
associated with a detectable level of genetic differentia- 
tion at nuclear markers. Such cryptic species have been 
identified in a diversity of organisms such as bryophytes, 
fungi, elasmobranches and arthropods [1,4,27,29-32]. 
Within Heliconius, a case of cryptic speciation has been 
recently discovered in two closely related species: H. 
timareta and H melpomene [33] . However, in H e. ches- 
tertonii we have found that despite deep divergence in 
mtDNA between two principal clades (2.15%), the AFLP 
analysis shows H. e. chestertonii to be a genetically 
homogeneous species (Figure 2). 

Endosymbionts infection 

Discordance between maternally inherited genetic mar- 
kers and those transmitted biparentally, can often be 
explained by the spread of endosymbionts such as Wol- 
bachia, Spiroplasma, Rickettsia, Arsenophonus, Cardi- 
nium, and others [34-37]. In most of the cases, these 
parasites are transmitted together with the mitochon- 
drial genome through the egg cytoplasm, so associations 
over time can be detected when the mitochondrial gen- 
ome is analyzed [38]. Some of these endosymbiont 
microorganisms can lead to reproductive alterations in 
their arthropod hosts and will lead to divergence 
between populations [19,39]. This is the first published 
study in which the presence of endosymbionts is tested 
in Heliconius, and we have identified two distinct 
lineages of Wolbachia within H e. chestertonii. However, 
neither of the principal mitochondrial clades of H. e. 
chestertonii is directly associated with the infection 
(Table 2). We did not find evidence of presence of other 
endosymbionts in populations of H. e. chestertonii. We 
should add a caveat to our results, which is that our 
PCR assay might not have detected all possible infected 
individuals. Indeed, the stage of development of the host 
can lead to over or underestimates of the real density of 
Wolbachia within populations [22]. This study was lim- 
ited to adult butterflies, and future analysis might 



include other stages such as egg, larvae and pupae. 
Nonetheless, the density of Wolbachia estimated within 
H e. chestertonii populations (7%) would be considered 
"very low" according to a recent classification [40]. In 
the future it would be interesting to further investigate 
the phenotypic effects of this infection. It is also of 
course possible that another endosymbiont is present, 
which was not sampled with the PCR assays described 
here. However, for the moment there is no evidence 
that the mitochondrial structure in H e. chestertonii is a 
result of endosymbiont infection. In summary, we have 
provided no evidence that mtDNA structure in H. e. 
chestertonii is due to either cryptic speciation or endo- 
symbiotic bacteria. This leaves historical processes 
within the species as the most likely cause for the 
pattern. 

Historical processes 

Climatic changes during glacial and interglacial periods 
can lead to contractions, expansions and fragmentations 
of populations [41-44]. The 2% divergence between 
lineages within H e. chestertonii suggests divergence 
within approximately the last million years. Even with- 
out such vicariance, isolation by distance can lead to 
genetic structure within species. Intriguingly, recent 
work on H cydno has shown that this species similarly 
has a marked genetic break in the center of its Cauca 
Valley range (Arias and Salazar, pers. comm.). In both 
species, southern Cauca populations are more closely 
related to subspecies on the Pacific coast than those in 
the North. A plausible scenario is that the Cauca Valley 
has been subject to a double colonization first from the 
central Andean valleys and subsequently from the Paci- 
fic populations in the west. Nonetheless, we are not 
aware of any geological evidence that would support 
this hypothesis. 

Conclusions 

Our genetic analysis shows that northern and southern 
populations of H e. chestertonii are genetically differen- 
tiated, but with only frequency differences in mtDNA 
clades, and no corresponding genetic structure at 
nuclear markers, despite the deep divergence (2.15%) 
between the two principal mtDNA clades. Our results 
support the assertion that mtDNA evidence alone 
should be used with caution in delimiting species 
boundaries. In this case, divergent haplotype groups 
within populations could not be explained by either 
cryptic speciation or endosymbiont infections. 

Methods 

Butterflies collections 

Between 1997 and 2008, 315 specimens of Heliconius 
erato chestertonii and H e. venus were collected at 15 
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different localities throughout southwest Colombia and 
one locality in Panama (Figure 1, Table 3). These two 
members of H. erato form a bimodal hybrid zone man- 
tained by strong premating isolation in one of these 
Colombian localities [17,26]. The body of each butterfly 
was separated from the wings and preserved in DMSO 
96% at -80°C, while wings were kept in glassine envel- 
opes. Tissue and wings are stored in the Instituto de 
Genetica, Universidad de los Andes, Colombia or in 
CDJs laboratory, University of Cambridge, UK (Table 3). 

DNA Butterflies extraction 

Genomic DNA was extracted from one-third of thorax 
or the end of the abdomen using DNeasy Blood & Tis- 
sue kit (Qiagen) following the manufacturer's tissue- 
extraction protocol. The thorax extractions were used 
for mitochondrial and AFLPs analysis and abdomen 
extractions for endosymbiont assays (Table 3). 

Sequencing and phylogenetic analysis of mitochondrial 
markers 

At least six individuals of each population, three of each 
sex when possible (Table 3), were used to amplify a 
total of 1551 bp of a mitochondrial region, which cov- 
ered the subunits I and II of Cytocrome Oxidase (CO 
I_II) and leucine-tRNA. We used published primers and 
PCR conditions from Beltran et al. 2002 for our amplifi- 
cations. Subsequently, all the PCR products were puri- 
fied for sequencing with Exonuclease I and Shrimp 
Alkaline Phosphatase enzymes (Fermentas) and sent to 



Macrogen Sequencing Service (Macrogen, Korea). 
Sequence editing was performed using Geneious Ver. 
5.4 [45]. Sequences were aligned using Clustal W Ver. 
2.0 [46] and checked for reading-frame errors in pro- 
tein-coding regions with MacClade Ver. 4.08 [47]. 
Sequences were deposited in GenBank (Access numbers: 
JF912810-JF912880). To complement our analysis, we 
included publicly available sequences of H. erato 
[13,16,17] and sequences of H. hecalesia and H. clysoni- 
mus were used as outgroups (GenBank access numbers 
in Additional file 5). Phylogenetic analyses were per- 
formed using two different approaches: Maximum-Parsi- 
mony (MP) and Bayesian inference (BI). 

Maximum-parsimony (MP) methods were imple- 
mented in PAUP 4.0b8 [48], using a heuristic search 
and tree bisection-reconnection (TBR) branch swap- 
ping option. A majority rule consensus tree was com- 
puted whenever multiple equally parsimonious trees 
were obtained. Parsimony bootstrap support values 
were estimated through 1000 bootstrap replicates. A 
Bayesian analysis was conducted with GTR + I + G 
nucleotide substitution model, which was the best-fit 
model obtained in JModel Test Ver. 0.1.1 [49] based 
on a Hierarchical Likelihood Ratio Test (HLRT). Baye- 
sian inference (BI) was carried out using four simulta- 
neous chains for ten million generations of Markov 
Chain Montecarlo (MCMC), sampling every 100 gen- 
erations. The consensus tree and posterior probability 
of the nodes was estimated with Mr Bayes Ver. 3.1 
[50,51]. 



Table 3 Populations collected and individuals used in each analysis 

Species Locality (Code) Latitude Longitud 





North 


West 


mtDNA 


AFLPs 


Endosymbionts 


H. e. chestertonii Marsella (MR) 


4° 52.39' 


75° 42.37' 


5/1 ¥ 


12 


33 


Trujillo CTR) 


4° 13.30' 


76° 20.30' 


7/2 ¥ 


11 


13 


Yotoco (YO) 


3° 51.12' 


76° 34.22' 


4 


4/6* 


4 


Buenos Aires (BA) 


3° 51.09' 


76° 25.37' 


6 


11 


18 


Carbonero (CR) 


3° 44.53' 


76° 28.56' 


6 


11 


14 


Caimital (CA) 


3° 44.12' 


76° 30.41' 


6 


11 


15 


Atuncela (AT) 


3° 44.05' 


76° 41.81' 


6 


12 


20 


Miravalle (MV) 


3° 41.22' 


76° 21.18' 


6 


11 


21 


Montahitas (MO) 


3° 41.03' 


76° 31.33' 


1/5 ¥ 


11 


20 


La Cumbre (CU) 


3° 38.04' 


76° 33.30' 


6 


11 


14 


Saladito (SA) 


3° 22.10' 


76° 39.04' 


4/2 ¥ 


9 


10 


Pance (PA) 


3° 19.27' 


76° 38.11' 


6 


10 


28 


Villa Colombia (VC) 


3° 11.51' 


76° 42.46' 


6 


11 


19 


H. e. chestertonii/CaWma River Valley 


3° 53.60' 


76° 37.57' 


2_5 ¥ /5 ¥ /2 ¥ 


11/11/2 


27/16 


H. e. venus /hybrids (CV) 












H. e. venus Juanchaco (JU) 


3° 56.22' 


77° 22.08' 


5 ¥ 


12 


29 


Rio Piedras (RP) + 


7° 63.62' 


78° 18.97' 




8* 





Sequences available and deposited by Arias et al. 2008 [GenBank: EU707581- EU707607]; individuals of Butterfly genetics group collection (Chris Jiggins 
Laboratory, Cambridge UK); + Darien Panama 
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Topology test and population genetics analysis 

The Shimodaira-Hasegawa (SH) log-likelihood test, as 
implemented in PAUP 4.0b8 [48,52], was used to test 
the monophyly of H. e. chestertonii. To test this a priori 
hypothesis, BI analysis was first performed using the 
same parameters as described earlier but this time by 
enforcing the monophyly of the complex as a topologi- 
cal constraint. The SH test was then used to compare 
trees obtained from both constrained and unconstrained 
analysis. 

We described DNA polymorphism in populations of 
H. e. chestertonii and H. e. venus (Table 1) with DnaSP 
Ver. 5.10.01 [53]. The measures employed were: 1) 
number of segregating sites («S), 2) number of haplotypes 
(h), 3) average pairwise number of differences between 
sequences (k) and genetic diversity (6). An analysis of 
molecular variance (AMOVA) and F S t values were 
obtained in ARLEQUIN Ver. 3.5 to determine genetic 
structure between populations from mtDNA markers 
[54]. To investigate relationships between populations a 
haplotype network was constructed using NETWORK 
Ver. 4.5.1.6 [55]. 

AFLPs, cluster analysis and genetic distance 

Genomic DNA from 185 individuals from all the popu- 
lations, was used for DNA fingerprinting with AFLP 
markers (Table 3). The DNA quality was quantified by 
spectrophotometry and only samples with an A260/ 
A280 ratio between 1.8 and 2.0 were used. For AFLP 
generation, we applied the general method of Vos et al. 
[56] with minor modifications. The AFLP® Core 
Reagent Kit (Invitrogen) was used for the digestion of 
125 ng of DNA for samples with EcoRl and Msel 
restriction endonucleases and the ligation with EcoRl/ 
Msel adaptors, following the manufacturers protocol. 
The best primer combinations (high variation and frag- 
ments between 50 to 400 bp) for pre- and selective 
amplifications were selected by comparing the final 
sequencing of eight different primer mixtures for 48 
individuals (Additional file 6). Sixteen samples chosen at 
random were run in duplicate during the screening pro- 
cess to ensure reproducibility of selected markers. As 
result, primers carrying a single selective nucleotide at 
the 3-end: £coRI+A (5-GAC TGC GTA CCA ATT CA- 
3) and Msel+C (5-GAT GAG TCC TGA GTA AC-3) 
were chosen for the pre-amplification cycle. The condi- 
tions in this first PCR were: denaturation at 94°C for 
60s, followed by 20 cycles with denaturation at 94°C for 
30s, annealing at 56°C for 60s and extension at 72°C for 
120s. Subsequently, selective amplifications were per- 
formed using four pairs of primers that contained three 
selective bases and EcoRl primers were labelled with the 
fluorescent dyes: EcoRl (FAM)+ACA/MseI+CGT, EcoRl 
(VIQ+ACC/Msel+CGT, EcoRl (PET)+ACT/MseI+CGT 



and EcoRl (NED)+ACG/MseI+CAC. For these amplifica- 
tions we applied a touchdown PCR of 12 cycles with 
denaturation at 94°C for 30s, an annealing temperature 
with 0.7°C stepwise reduction from 65°C to 56°C for 60s 
and extension at 72°C for 120s, after which 23 additional 
cycles were run with fixed annealing temperature of 56° 
C for 30s. AFLP reaction products were labelled with an 
internal size standard (GeneScan™-500 LIZ™) and run 
on an ABI 3730 automated DNA sequencer (Applied 
Biosystems). Fragment data were collected and analysed 
with GeneMapper Ver. 4.0 (Applied Biosystems). All 
bands were visually confirmed and only unambiguous 
loci between 50 and 400 bp were included in the analy- 
sis. Each individual genotype was manually assessed. A 
binary data matrix of presence (1) or absence (0) of 
bands of each size by sample was generated. 

We used the presence/absence matrix of AFLP frag- 
ments to differentiate all individuals and populations 
using a Bayesian clustering algorithm implemented in 
the program Structure Ver. 2.3.3 [57,58]. The number of 
clusters {K) was determined by comparing the likelihood 
ratios for K values between 1 and 6. Previous runs 
including more than 6 clusters were done and the likeli- 
hood values were lower and are not shown in this analy- 
sis. Each likelihood value was estimated with runs that 
involved a burn-in of between 10 4 and 10 7 MCMC gen- 
erations with 10 interactions. An admixture model 
under the assumption of Hardy- Weinberg equilibrium 
was implemented in the runs. The best number of clus- 
ters was confirmed with the ad hoc statistic AK [59]. 
Additionally we implemented a principal component 
analysis (PCA) using the software Genetix Ver. 4.05. 
[60] to confirm the results and variation of our data. 
Measures of heterozygosity and variation were estimated 
with AFLP-SURV Ver.1.0 [61] and with an AMOVA in 
Arlequin [54]. 

Determination of the presence of endosymbionts 

A total of 301 individuals (at least 10 from each popula- 
tion), were used to search for the presence of endosym- 
bionts (Table 2 and 3). We used PCR to test for the 
presence of Wolbachia, Ricketsia and Spiroplasma in H. 
e. chestertonii and H. e. venus. Where the presence of 
Wolbachia was detected, additional PCR was carried out 
to determine the Wolbachia clade (A or B). Amplifica- 
tions were carried out with specific primers that ampli- 
fied the wsp gene of Wolbachia (A and B-clade) and 
16S rDNA for the other two bacteria (Table 4). We 
designed three novel specific primers from our Helico- 
nius wsp gene sequences in addition to those proposed 
by Zhou et al. 1998 (Table 4). Different combinations of 
these primers allowed us to determine the presence of 
Wolbachia (diagnostic combination: H81F/H691R), Wol- 
bachia A-clade (H308F/H691R) and Wolbachia B-clade 
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Table 4 Primers for mtDNA and endosymbionts genes 


Gene 


Target species 


Primers 


Sequence {5'-3') 


Reference 


COI 


Heliconius 


Jerry 


CAACA I I A I I I I GA I I I I I I GG 


Beltran et al. 2002 






Pat 


^CAATGCACTAATCTGCCATA^A 




tRNA-leu 


Heliconius 


George 


TAGGATOGCTGGAATACC 


Beltran et al. 2002 


and COII 




Imelda 


CATOGAAGTAATOCTAAmAACTA 




Wsp 


Wolbachia 


H81F 


AACTAGCTACTACGTOGmGC 


This study 






H681R 


GCTACTCCAGOTCTGCAC 






A-clade 


H308R 


TOAAGATGTAACAmGACC 






B-clade 


wsp522R 


ACCAGI I I I IGOTGATA 


Zhou et al. 1998 


16S 


Ricketsia 


RSSUF 


CGGCmCAAAACTACTAATCTA 


von der Schulenburg 






RSSUR 


GAAAGCATCTCTGCGATCCG 


et al. 2001 


16S 


Spiroplasmo 


27f 


GAGAGmGATCCTGGCTCAG 


Hurst et al. 1999 






MGSO 


TGCACCATCTGTCACTCTGTOACCTC 


Van kuppeveld et a/. 1992 



(H81F/wsp522R). Positive controls for these PCR reac- 
tions were DNA from known infected insects supplied 
by Emily Hornett {Wolbachia) and Francis Jiggins (Rick- 
etsia and Spiroplasmo) in the University of Cambridge, 
UK. We followed the amplification conditions as 
described previously [35,62-65]. 

Additional material 
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